      PROGRAM adv1d

      use csv_file; !!

      implicit none

!      INTEGER :: i0=100
      INTEGER :: i0
      INTEGER :: n,i
      parameter (n=400)
      DOUBLE PRECISION, DIMENSION(n) :: c(n)
      DOUBLE PRECISION, DIMENSION(n) :: x(n)
      DOUBLE PRECISION, DIMENSION(n) :: z,z1,z1_d,z2,z2_d
      data z /n*1./
      data z1 /n*1.0/
      data z2 /n*1.1/

!    Open our output file
      call system('rm plot1d.csv')
      open(16, file = 'plot1d.csv') !!

!    Get i0 (no idea why)
      print *,'enter i0'
!      read *,i0
       i0 = 5
      print *,'c=',360.0/float(i0)

!    Loop through the grid?
      do 1000 i=1,n
          c(i)=360./float(i0)
 1000     x(i)=float(i)
          z(372)=5.
          z(373)=5.
          z(374)=5.
!         write (16) x,z
          z1(12)=3.0
          z1(13)=3.0
          z1(14)=3.0
          z2(12)=3.3
          z2(13)=3.3
          z2(14)=3.3

!         Writing original data
          call csv_write(16, 'z1_orig', .false.) !!
          call csv_write(16, z1, .true.) !! 
!         Writing original data
          call csv_write(16, 'z2_orig', .false.) !!
          call csv_write(16, z2, .true.) !!

          do 1020 i=1,i0
!              call adv4p (z1,c)
!              call adv4m (z1,c)
              call adv4m_d (z1, z1_d, c)
              call adv4m_d (z2, z2_d, c)
 1020     continue
!     print *, z1
!      write (16) z1,z2

!   Writing to a csv instead
!          PRINT '(F8.3,A2,F8.3)', z1,', ',z2
!         WRITE (16, '(F8.3, F8.3)'), z1, z2

!			Writing z1
          call csv_write(16, 'z1', .false.) !!
          call csv_write(16, z1, .true.) !!
          call csv_write(16, 'z1_d', .false.) !!
          call csv_write(16, z1_d, .true.) !!

!			Writing z2
          call csv_write(16, 'z2', .false.) !!
          call csv_write(16, z2, .true.) !!
          call csv_write(16, 'zd_2', .false.) !!
          call csv_write(16, z2_d, .true.) !!

           close(16) !!
      stop
      end

      subroutine adv0 (y,c)
      implicit none
! upstream procedure

      INTEGER :: n,i
      parameter (n=400)
      DOUBLE PRECISION, DIMENSION(n) :: fm,fp,y,c
      do 1000 i=1,n-1
          fm(i)=-amin1(0.,c(i))*y(i+1)
          fp(i)=amax1(0.,c(i))*y(i)
 1000 continue
      do 1010 i=2,n-1
          y(i)=y(i)-fm(i-1)+fp(i-1)+fm(i)-fp(i)
!          print  *, "y(",i,") = ", y(i-1)
 1010 continue
      print *, " "
      return

      end

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.1 (r2754) - 01/12/2009 09:44
!  
!  Differentiation of adv0 in forward (tangent) mode:
!   variations  of output variables: y
!   with respect to input variables: y
      SUBROUTINE ADV0_D0(y, yd, c)
        IMPLICIT NONE
      ! upstream procedure
        INTEGER :: n, i
        PARAMETER (n=400)
        DOUBLE PRECISION, DIMENSION(n) :: fm, fp, y, c
        DOUBLE PRECISION, DIMENSION(n) :: fmd, fpd, yd
        DOUBLE PRECISION :: amin11
        DOUBLE PRECISION :: amax11
        INTRINSIC AMAX1
        INTRINSIC AMIN1
        fmd(1:n) = 0.0
        fpd(1:n) = 0.0
        DO i=1,n-1
          IF (0. .GT. c(i)) THEN
            amin11 = c(i)
          ELSE
            amin11 = 0.
          END IF
          fmd(i) = -(amin11*yd(i+1))
          fm(i) = -(amin11*y(i+1))
          IF (0. .LT. c(i)) THEN
            amax11 = c(i)
          ELSE
            amax11 = 0.
          END IF
          fpd(i) = amax11*yd(i)
          fp(i) = amax11*y(i)
        END DO
        DO i=2,n-1
          yd(i) = yd(i) - fmd(i-1) + fpd(i-1) + fmd(i) - fpd(i)
          y(i) = y(i) - fm(i-1) + fp(i-1) + fm(i) - fp(i)
        END DO
        RETURN
      END SUBROUTINE ADV0_D0


      subroutine adv2p (y,c)
! area preserving flux form; Bott (1989): Monthly Weather Review.
! second order positive definite version.
! y(i) is transport quantity, input and output.
! boundary conditions are y(1)=const, y(n)=const.
! c(i) is Courant number satisfying the CFL criterion, input.
! fm(i), fp(i) are fluxes for u(i)<0 and u(i)>0, respectively.
! a0, a1, a2, a3, a4 are coefficients of polynomials in gridbox i.
! At i=1 and i=n first order polynomial,
! at 2<=i<=n-1 second order polynomial.
! w(i) are weighting factors.
! the numerical grid is equidistant.
! the procedure is one dimensional.
! for multidimensional applications time splitting has to be used.
! the quantities c(i), fm(i), fp(i)  are given at the right
! boundary of grid cell i.
! Thus, fm(i) is flux from gridbox i+1 into gridbox i for c(i)<0,
! fp(i) is flux from gridbox i into gridbox i+1 for c(i)>0.
      parameter (n=400)
      dimension fm(n),fp(n),w(n),y(n),c(n)
      cr=amax1(0.,c(1))
      fp(1)=amin1(y(1),cr*(y(1)+(1.-cr)*(y(2)-y(1))*0.5))
      w(1)=1.
      do 1000 i=2,n-1
      a0=(26.*y(i)-y(i+1)-y(i-1))/24.
      a1=(y(i+1)-y(i-1))/16.
      a2=(y(i+1)+y(i-1)-2.*y(i))/48.
      cl=-amin1(0.,c(i-1))
      x1=1.-2.*cl
      x2=x1*x1
      fm(i-1)=amax1(0.,a0*cl-a1*(1.-x2)+a2*(1.-x1*x2))
      cr=amax1(0.,c(i))
      x1=1.-2.*cr
      x2=x1*x1
      fp(i)=amax1(0.,a0*cr+a1*(1.-x2)+a2*(1.-x1*x2))
      w(i)=y(i)/amax1(fm(i-1)+fp(i)+1.e-15,a0+2.*a2)
 1000 continue
      cl=-amin1(0.,c(n-1))
      fm(n-1)=amin1(y(n),cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5))
      w(n)=1.
      do 1010 i=2,n-1
 1010 y(i)=y(i)-(fm(i-1)+fp(i))*w(i)+fm(i)*w(i+1)+fp(i-1)*w(i-1)
      return
      end

      subroutine adv4p (y,c)
! area preserving flux form; Bott (1989): Monthly Weather Review.
! fourth order positive definite version.
! y(i) is transport quantity, input and output.
! boundary conditions are y(1)=const, y(n)=const.
! c(i) is Courant number satisfying the CFL criterion, input.
! fm(i), fp(i) are fluxes for u(i)<0 and u(i)>0, respectively.
! a0, a1, a2, a3, a4 are coefficients of polynomials in gridbox i.
! At i=1 and i=n first order polynomial,
! at i=2 and i=n-1 second order polynomial,
! at 3<=i<=n-2 fourth order polynomial.
! w(i) are weighting factors.
! the numerical grid is equidistant.
! the procedure is one dimensional.
! for multidimensional applications time splitting has to be used.
! the quantities c(i), fm(i), fp(i)  are given at the right
! boundary of grid cell i.
! Thus, fm(i) is flux from gridbox i+1 into gridbox i for c(i)<0,
! fp(i) is flux from gridbox i into gridbox i+1 for c(i)>0.

      INTEGER :: n

      parameter (n=400)
      dimension y(n),c(n),fm(n),fp(n),w(n)
      cr=amax1(0.,c(1))
      fp(1)=amin1(y(1),cr*(y(1)+(1.-cr)*(y(2)-y(1))*0.5))
      w(1)=1.
      a0=(26.*y(2)-y(3)-y(1))/24.
      a1=(y(3)-y(1))/16.
      a2=(y(3)+y(1)-2.*y(2))/48.
      cl=-amin1(0.,c(1))
      x1=1.-2.*cl
      x2=x1*x1
      fm(1)=amax1(0.,a0*cl-a1*(1.-x2)+a2*(1.-x1*x2))
      cr=amax1(0.,c(2))
      x1=1.-2.*cr
      x2=x1*x1
      fp(2)=amax1(0.,a0*cr+a1*(1.-x2)+a2*(1.-x1*x2))
      w(2)=y(2)/amax1(fm(1)+fp(2)+1.e-15,a0+2.*a2)
      do 1000 i=3,n-2
          a0=(9.*(y(i+2)+y(i-2))-116.*(y(i+1)+y(i-1))+2134.*y(i))/1920.0
          a1=(-5.*(y(i+2)-y(i-2))+34.*(y(i+1)-y(i-1)))/384.
          a2=(-y(i+2)+12.*(y(i+1)+y(i-1))-22.*y(i)-y(i-2))/384.
          a3=(y(i+2)-2.*(y(i+1)-y(i-1))-y(i-2))/768.
          a4=(y(i+2)-4.*(y(i+1)+y(i-1))+6.*y(i)+y(i-2))/3840.
          cl=-amin1(0.,c(i-1))
          x1=1.-2.*cl
          x2=x1*x1
          x3=x1*x2
          fm(i-1)=amax1(0.,a0*cl-a1*(1.-x2)+a2*(1.-x3)-a3*(1.-x1*x3)+a4*(1.-x2*x3))
          cr=amax1(0.,c(i))
          x1=1.-2.*cr
          x2=x1*x1
          x3=x1*x2
          fp(i)=amax1(0.,a0*cr+a1*(1.-x2)+a2*(1.-x3)+a3*(1.-x1*x3)+a4*(1.-x2*x3))
          w(i)=y(i)/amax1(fm(i-1)+fp(i)+1.e-15,y(i))
 1000 continue
      a0=(26.*y(n-1)-y(n)-y(n-2))/24.
      a1=(y(n)-y(n-2))/16.
      a2=(y(n)+y(n-2)-2.*y(n-1))/48.
      cl=-amin1(0.,c(n-2))
      x1=1.-2.*cl
      x2=x1*x1
      fm(n-2)=amax1(0.,a0*cl-a1*(1.-x2)+a2*(1.-x1*x2))
      cr=amax1(0.,c(n-1))
      x1=1.-2.*cr
      x2=x1*x1
      fp(n-1)=amax1(0.,a0*cr+a1*(1.-x2)+a2*(1.-x1*x2))
      w(n-1)=y(n-1)/amax1(fm(n-2)+fp(n-1)+1.e-15,a0+2.*a2)
      cl=-amin1(0.,c(n-1))
      fm(n-1)=amin1(y(n),cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5))
      w(n)=1.
      do 1010 i=2,n-1
 1010 y(i)=y(i)-(fm(i-1)+fp(i))*w(i)+fm(i)*w(i+1)+fp(i-1)*w(i-1)
      return
      end

      !        Generated by TAPENADE     (INRIA, Tropics team)
      !  Tapenade 3.1 (r2754) - 01/12/2009 09:44
      !  
      !  Differentiation of adv4p in forward (tangent) mode:
      !   variations  of output variables: y
      !   with respect to input variables: y
      SUBROUTINE ADV4P_D(y, yd, c)
        IMPLICIT NONE
      ! area preserving flux form; Bott (1989): Monthly Weather Review.
      ! fourth order positive definite version.
      ! y(i) is transport quantity, input and output.
      ! boundary conditions are y(1)=const, y(n)=const.
      ! c(i) is Courant number satisfying the CFL criterion, input.
      ! fm(i), fp(i) are fluxes for u(i)<0 and u(i)>0, respectively.
      ! a0, a1, a2, a3, a4 are coefficients of polynomials in gridbox i.
      ! At i=1 and i=n first order polynomial,
      ! at i=2 and i=n-1 second order polynomial,
      ! at 3<=i<=n-2 fourth order polynomial.
      ! w(i) are weighting factors.
      ! the numerical grid is equidistant.
      ! the procedure is one dimensional.
      ! for multidimensional applications time splitting has to be used.
      ! the quantities c(i), fm(i), fp(i)  are given at the right
      ! boundary of grid cell i.
      ! Thus, fm(i) is flux from gridbox i+1 into gridbox i for c(i)<0,
      ! fp(i) is flux from gridbox i into gridbox i+1 for c(i)>0.
        INTEGER :: n
        PARAMETER (n=400)
        DOUBLE PRECISION :: w
        DOUBLE PRECISION :: wd
        DOUBLE PRECISION :: fm
        DOUBLE PRECISION :: fmd
        DOUBLE PRECISION :: fp
        DOUBLE PRECISION :: fpd
        DOUBLE PRECISION :: y
        DOUBLE PRECISION :: c
        DIMENSION y(n), c(n), fm(n), fp(n), w(n)
        DOUBLE PRECISION :: yd
        DIMENSION yd(n), fmd(n), fpd(n), wd(n)
        DOUBLE PRECISION :: cr
        DOUBLE PRECISION :: a0
        DOUBLE PRECISION :: a0d
        DOUBLE PRECISION :: a1
        DOUBLE PRECISION :: a1d
        DOUBLE PRECISION :: a2
        DOUBLE PRECISION :: a2d
        DOUBLE PRECISION :: cl
        DOUBLE PRECISION :: x1
        DOUBLE PRECISION :: x2
        INTEGER :: i
        DOUBLE PRECISION :: a3
        DOUBLE PRECISION :: a3d
        DOUBLE PRECISION :: a4
        DOUBLE PRECISION :: a4d
        DOUBLE PRECISION :: x3
        DOUBLE PRECISION :: amax12d
        DOUBLE PRECISION :: y4d
        DOUBLE PRECISION :: amin14
        DOUBLE PRECISION :: amin13
        DOUBLE PRECISION :: amin12
        DOUBLE PRECISION :: amin11
        DOUBLE PRECISION :: amax11d
        DOUBLE PRECISION :: y3d
        DOUBLE PRECISION :: amax13
        DOUBLE PRECISION :: amax12
        DOUBLE PRECISION :: amax11
        INTRINSIC AMAX1
        DOUBLE PRECISION :: y2d
        DOUBLE PRECISION :: amax13d
        INTRINSIC AMIN1
        DOUBLE PRECISION :: y4
        DOUBLE PRECISION :: y3
        DOUBLE PRECISION :: y2
        DOUBLE PRECISION :: y1
        DOUBLE PRECISION :: y1d
        IF (0. .LT. c(1)) THEN
          cr = c(1)
        ELSE
          cr = 0.
        END IF
        y1d = cr*(yd(1)+(1.-cr)*0.5*(yd(2)-yd(1)))
        y1 = cr*(y(1)+(1.-cr)*(y(2)-y(1))*0.5)
        IF (y(1) .GT. y1) THEN
          fpd(1) = y1d
          fp(1) = y1
        ELSE
          fpd(1) = yd(1)
          fp(1) = y(1)
        END IF
        wd(1) = 0.0
        w(1) = 1.
        a0d = (26.*yd(2)-yd(3)-yd(1))/24.
        a0 = (26.*y(2)-y(3)-y(1))/24.
        a1d = (yd(3)-yd(1))/16.
        a1 = (y(3)-y(1))/16.
        a2d = (yd(3)+yd(1)-2.*yd(2))/48.
        a2 = (y(3)+y(1)-2.*y(2))/48.
        IF (0. .GT. c(1)) THEN
          amin11 = c(1)
        ELSE
          amin11 = 0.
        END IF
        cl = -amin11
        x1 = 1. - 2.*cl
        x2 = x1*x1
        IF (0. .LT. a0*cl - a1*(1.-x2) + a2*(1.-x1*x2)) THEN
          fmd(1) = cl*a0d - (1.-x2)*a1d + (1.-x1*x2)*a2d
          fm(1) = a0*cl - a1*(1.-x2) + a2*(1.-x1*x2)
        ELSE
          fmd(1) = 0.0
          fm(1) = 0.
          fmd(1:n) = 0.0
        END IF
        IF (0. .LT. c(2)) THEN
          cr = c(2)
        ELSE
          cr = 0.
        END IF
        x1 = 1. - 2.*cr
        x2 = x1*x1
        y2d = cr*a0d + (1.-x2)*a1d + (1.-x1*x2)*a2d
        y2 = a0*cr + a1*(1.-x2) + a2*(1.-x1*x2)
        IF (0. .LT. y2) THEN
          fpd(2) = y2d
          fp(2) = y2
        ELSE
          fpd(2) = 0.0
          fp(2) = 0.
        END IF
        IF (fm(1) + fp(2) + 1.e-15 .LT. a0 + 2.*a2) THEN
          amax11d = a0d + 2.*a2d
          amax11 = a0 + 2.*a2
        ELSE
          amax11d = fmd(1) + fpd(2)
          amax11 = fm(1) + fp(2) + 1.e-15
        END IF
        wd(2) = (yd(2)*amax11-y(2)*amax11d)/amax11**2
        w(2) = y(2)/amax11
        DO i=3,n-2
          a0d = (9.*(yd(i+2)+yd(i-2))-116.*(yd(i+1)+yd(i-1))+2134.*yd(i))/&
      &      1920.0
          a0 = (9.*(y(i+2)+y(i-2))-116.*(y(i+1)+y(i-1))+2134.*y(i))/1920.0
          a1d = (34.*(yd(i+1)-yd(i-1))-5.*(yd(i+2)-yd(i-2)))/384.
          a1 = (-(5.*(y(i+2)-y(i-2)))+34.*(y(i+1)-y(i-1)))/384.
          a2d = (12.*(yd(i+1)+yd(i-1))-yd(i+2)-22.*yd(i)-yd(i-2))/384.
          a2 = (-y(i+2)+12.*(y(i+1)+y(i-1))-22.*y(i)-y(i-2))/384.
          a3d = (yd(i+2)-2.*(yd(i+1)-yd(i-1))-yd(i-2))/768.
          a3 = (y(i+2)-2.*(y(i+1)-y(i-1))-y(i-2))/768.
          a4d = (yd(i+2)-4.*(yd(i+1)+yd(i-1))+6.*yd(i)+yd(i-2))/3840.
          a4 = (y(i+2)-4.*(y(i+1)+y(i-1))+6.*y(i)+y(i-2))/3840.
          IF (0. .GT. c(i-1)) THEN
            amin12 = c(i-1)
          ELSE
            amin12 = 0.
          END IF
          cl = -amin12
          x1 = 1. - 2.*cl
          x2 = x1*x1
          x3 = x1*x2
          IF (0. .LT. a0*cl - a1*(1.-x2) + a2*(1.-x3) - a3*(1.-x1*x3) + a4*(1.&
      &        -x2*x3)) THEN
            fmd(i-1) = cl*a0d - (1.-x2)*a1d + (1.-x3)*a2d - (1.-x1*x3)*a3d + (&
      &        1.-x2*x3)*a4d
            fm(i-1) = a0*cl - a1*(1.-x2) + a2*(1.-x3) - a3*(1.-x1*x3) + a4*(1.&
      &        -x2*x3)
          ELSE
            fmd(i-1) = 0.0
            fm(i-1) = 0.
          END IF
          IF (0. .LT. c(i)) THEN
            cr = c(i)
          ELSE
            cr = 0.
          END IF
          x1 = 1. - 2.*cr
          x2 = x1*x1
          x3 = x1*x2
          y3d = cr*a0d + (1.-x2)*a1d + (1.-x3)*a2d + (1.-x1*x3)*a3d + (1.-x2*&
      &      x3)*a4d
          y3 = a0*cr + a1*(1.-x2) + a2*(1.-x3) + a3*(1.-x1*x3) + a4*(1.-x2*x3)
          IF (0. .LT. y3) THEN
            fpd(i) = y3d
            fp(i) = y3
          ELSE
            fpd(i) = 0.0
            fp(i) = 0.
          END IF
          IF (fm(i-1) + fp(i) + 1.e-15 .LT. y(i)) THEN
            amax12d = yd(i)
            amax12 = y(i)
          ELSE
            amax12d = fmd(i-1) + fpd(i)
            amax12 = fm(i-1) + fp(i) + 1.e-15
          END IF
          wd(i) = (yd(i)*amax12-y(i)*amax12d)/amax12**2
          w(i) = y(i)/amax12
        END DO
        a0d = (26.*yd(n-1)-yd(n)-yd(n-2))/24.
        a0 = (26.*y(n-1)-y(n)-y(n-2))/24.
        a1d = (yd(n)-yd(n-2))/16.
        a1 = (y(n)-y(n-2))/16.
        a2d = (yd(n)+yd(n-2)-2.*yd(n-1))/48.
        a2 = (y(n)+y(n-2)-2.*y(n-1))/48.
        IF (0. .GT. c(n-2)) THEN
          amin13 = c(n-2)
        ELSE
          amin13 = 0.
        END IF
        cl = -amin13
        x1 = 1. - 2.*cl
        x2 = x1*x1
        IF (0. .LT. a0*cl - a1*(1.-x2) + a2*(1.-x1*x2)) THEN
          fmd(n-2) = cl*a0d - (1.-x2)*a1d + (1.-x1*x2)*a2d
          fm(n-2) = a0*cl - a1*(1.-x2) + a2*(1.-x1*x2)
        ELSE
          fmd(n-2) = 0.0
          fm(n-2) = 0.
        END IF
        IF (0. .LT. c(n-1)) THEN
          cr = c(n-1)
        ELSE
          cr = 0.
        END IF
        x1 = 1. - 2.*cr
        x2 = x1*x1
        y4d = cr*a0d + (1.-x2)*a1d + (1.-x1*x2)*a2d
        y4 = a0*cr + a1*(1.-x2) + a2*(1.-x1*x2)
        IF (0. .LT. y4) THEN
          fpd(n-1) = y4d
          fp(n-1) = y4
        ELSE
          fpd(n-1) = 0.0
          fp(n-1) = 0.
        END IF
        IF (fm(n-2) + fp(n-1) + 1.e-15 .LT. a0 + 2.*a2) THEN
          amax13d = a0d + 2.*a2d
          amax13 = a0 + 2.*a2
        ELSE
          amax13d = fmd(n-2) + fpd(n-1)
          amax13 = fm(n-2) + fp(n-1) + 1.e-15
        END IF
        wd(n-1) = (yd(n-1)*amax13-y(n-1)*amax13d)/amax13**2
        w(n-1) = y(n-1)/amax13
        IF (0. .GT. c(n-1)) THEN
          amin14 = c(n-1)
        ELSE
          amin14 = 0.
        END IF
        cl = -amin14
        IF (y(n) .GT. cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5)) THEN
          fmd(n-1) = cl*(yd(n)-(1.-cl)*0.5*(yd(n)-yd(n-1)))
          fm(n-1) = cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5)
        ELSE
          fmd(n-1) = yd(n)
          fm(n-1) = y(n)
        END IF
        wd(n) = 0.0
        w(n) = 1.
        DO i=2,n-1
          yd(i) = yd(i) - (fmd(i-1)+fpd(i))*w(i) - (fm(i-1)+fp(i))*wd(i) + fmd&
      &      (i)*w(i+1) + fm(i)*wd(i+1) + fpd(i-1)*w(i-1) + fp(i-1)*wd(i-1)
          y(i) = y(i) - (fm(i-1)+fp(i))*w(i) + fm(i)*w(i+1) + fp(i-1)*w(i-1)
        END DO
        RETURN
      END SUBROUTINE ADV4P_D


      subroutine adv4m (y,c)
! area preserving flux form; Bott (1989): Monthly Weather Review.
! fourth order monotone version.
! y(i) is transport quantity, input and output.
! boundary conditions are y(1)=const, y(n)=const.
! c(i) is Courant number satisfying the CFL criterion, input.
! fm(i), fp(i) are fluxes for u(i)<0 and u(i)>0, respectively.
! a0, a1, a2, a3, a4 are coefficients of polynomials in gridbox i.
! At i=1 and i=n first order polynomial,
! at i=2 and i=n-1 second order polynomial,
! at 3<=i<=n-2 fourth order polynomial.
! w(i) are weighting factors.
! the numerical grid is equidistant.
! the procedure is one dimensional.
! for multidimensional applications time splitting has to be used.
! the quantities c(i), fm(i), fp(i)  are given at the right
! boundary of grid cell i.
! Thus, fm(i) is flux from gridbox i+1 into gridbox i for c(i)<0,
! fp(i) is flux from gridbox i into gridbox i+1 for c(i)>0.

      INTEGER :: n

      parameter (n=400)
      DOUBLE PRECISION, DIMENSION(n) :: a0,a1,a2,a3,a4,y,c,fm,fp,w
      a0(2)=(26.*y(2)-y(3)-y(1))/24.
      a1(2)=(y(3)-y(1))/16.
      a2(2)=(y(3)+y(1)-2.*y(2))/48.
      a3(2)=0.
      a4(2)=0.
      do 1000 i=3,n-2
          a0(i)=(9.*(y(i+2)+y(i-2))-116.*(y(i+1)+y(i-1))+2134.*y(i))/1920.
          a1(i)=(-5.*(y(i+2)-y(i-2))+34.*(y(i+1)-y(i-1)))/384.
          a2(i)=(-y(i+2)+12.*(y(i+1)+y(i-1))-22.*y(i)-y(i-2))/384.
          a3(i)=(y(i+2)-2.*(y(i+1)-y(i-1))-y(i-2))/768.
 1000 a4(i)=(y(i+2)-4.*(y(i+1)+y(i-1))+6.*y(i)+y(i-2))/3840.
      a0(n-1)=(26.*y(n-1)-y(n)-y(n-2))/24.
      a1(n-1)=(y(n)-y(n-2))/16.
      a2(n-1)=(y(n)+y(n-2)-2.*y(n-1))/48.
      a3(n-1)=0.
      a4(n-1)=0.
      w(1)=1.
      w(n)=1.
      fm(n)=0.
      cl=-amin1(0.,c(n-1))
      fm(n-1)=amin1(y(n),cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5))
      clm=cl
      do 1020 i=n-1,2,-1
          cl=clm
          clm=-amin1(0.,c(i-1))
          x1=1.-2.*cl
          x2=x1*x1
          x3=x1*x2
          ymin=amin1(y(i),y(i+1))
          ymax=amax1(y(i),y(i+1))
          fmim=amax1(0.,a0(i)*cl-a1(i)*(1.-x2)+a2(i)*(1.-x3)-a3(i)*(1.-x1*x3)+a4(i)*(1.-x2*x3))
          fmim=amin1(fmim,y(i)-ymin+fm(i))
          fmim=amax1(fmim,y(i)-ymax+fm(i))
          fm(i-1)=amax1(0.,fmim-(cl-clm)*y(i))
 1020 continue
      cr=amax1(0.,c(1))
      fp(1)=amin1(y(1),cr*(y(1)+(1.-cr)*(y(2)-y(1))*0.5))
      crp=cr
      do 1030 i=2,n-1
          cr=crp
          crp=amax1(0.,c(i))
          x1=1.-2.*cr
          x2=x1*x1
          x3=x1*x2
          ymin=amin1(y(i-1),y(i))
          ymax=amax1(y(i-1),y(i))
          fpi=amax1(0.,a0(i)*cr+a1(i)*(1.-x2)+a2(i)*(1.-x3)+a3(i)*(1.-x1*x3)+a4(i)*(1.-x2*x3))
          fpi=amin1(fpi,y(i)-ymin+fp(i-1))
          fpi=amax1(fpi,y(i)-ymax+fp(i-1))
          fp(i)=amax1(0.,fpi-(cr-crp)*y(i))
          w(i)=y(i)/amax1(fm(i-1)+fp(i)+1.e-15,y(i))
 1030 continue
      do 1040 i=2,n-1
 1040 y(i)=y(i)-(fm(i-1)+fp(i))*w(i)+fm(i)*w(i+1)+fp(i-1)*w(i-1)
      return
      end

      !        Generated by TAPENADE     (INRIA, Tropics team)
      !  Tapenade 3.1 (r2754) - 01/12/2009 09:44
      !  
      !  Differentiation of adv4m in forward (tangent) mode:
      !   variations  of output variables: y
      !   with respect to input variables: y
      SUBROUTINE ADV4M_D(y, yd, c)
        IMPLICIT NONE
      ! area preserving flux form; Bott (1989): Monthly Weather Review.
      ! fourth order monotone version.
      ! y(i) is transport quantity, input and output.
      ! boundary conditions are y(1)=const, y(n)=const.
      ! c(i) is Courant number satisfying the CFL criterion, input.
      ! fm(i), fp(i) are fluxes for u(i)<0 and u(i)>0, respectively.
      ! a0, a1, a2, a3, a4 are coefficients of polynomials in gridbox i.
      ! At i=1 and i=n first order polynomial,
      ! at i=2 and i=n-1 second order polynomial,
      ! at 3<=i<=n-2 fourth order polynomial.
      ! w(i) are weighting factors.
      ! the numerical grid is equidistant.
      ! the procedure is one dimensional.
      ! for multidimensional applications time splitting has to be used.
      ! the quantities c(i), fm(i), fp(i)  are given at the right
      ! boundary of grid cell i.
      ! Thus, fm(i) is flux from gridbox i+1 into gridbox i for c(i)<0,
      ! fp(i) is flux from gridbox i into gridbox i+1 for c(i)>0.
        INTEGER :: n
        DOUBLE PRECISION :: c
        PARAMETER (n=400)
        DOUBLE PRECISION :: w
        DOUBLE PRECISION :: wd
        DOUBLE PRECISION :: a0
        DOUBLE PRECISION :: a0d
        DOUBLE PRECISION :: a1
        DOUBLE PRECISION :: a1d
        DOUBLE PRECISION :: a2
        DOUBLE PRECISION :: a2d
        DOUBLE PRECISION :: a3
        DOUBLE PRECISION :: a3d
        DOUBLE PRECISION :: a4
        DOUBLE PRECISION :: a4d
        DOUBLE PRECISION :: fm
        DOUBLE PRECISION :: fmd
        DOUBLE PRECISION :: fp
        DOUBLE PRECISION :: fpd
        DOUBLE PRECISION :: y
        DIMENSION a0(n), a1(n), a2(n), a3(n), a4(n), y(n), c(n), fm(n), fp(n)&
      &      , w(n)
        DOUBLE PRECISION :: yd
        DIMENSION a0d(n), a1d(n), a2d(n), a3d(n), a4d(n), yd(n), fmd(n), fpd(n&
      &      ), wd(n)
        INTEGER :: i
        DOUBLE PRECISION :: cl
        DOUBLE PRECISION :: clm
        DOUBLE PRECISION :: x1
        DOUBLE PRECISION :: x2
        DOUBLE PRECISION :: x3
        DOUBLE PRECISION :: ymin
        DOUBLE PRECISION :: ymind
        DOUBLE PRECISION :: ymax
        DOUBLE PRECISION :: ymaxd
        DOUBLE PRECISION :: fmim
        DOUBLE PRECISION :: fmimd
        DOUBLE PRECISION :: cr
        DOUBLE PRECISION :: crp
        DOUBLE PRECISION :: fpi
        DOUBLE PRECISION :: fpid
        DOUBLE PRECISION :: amin12
        DOUBLE PRECISION :: amin11
        DOUBLE PRECISION :: amax11d
        DOUBLE PRECISION :: amax11
        INTRINSIC AMAX1
        DOUBLE PRECISION :: y2d
        INTRINSIC AMIN1
        DOUBLE PRECISION :: y2
        DOUBLE PRECISION :: y1
        DOUBLE PRECISION :: y1d
        a0d(2) = (26.*yd(2)-yd(3)-yd(1))/24.
        a0(2) = (26.*y(2)-y(3)-y(1))/24.
        a1d(2) = (yd(3)-yd(1))/16.
        a1(2) = (y(3)-y(1))/16.
        a2d(2) = (yd(3)+yd(1)-2.*yd(2))/48.
        a2(2) = (y(3)+y(1)-2.*y(2))/48.
        a3d(2) = 0.0
        a3(2) = 0.
        a4d(2) = 0.0
        a4(2) = 0.
        a3d(1:n) = 0.0
        a4d(1:n) = 0.0
        DO i=3,n-2
          a0d(i) = (9.*(yd(i+2)+yd(i-2))-116.*(yd(i+1)+yd(i-1))+2134.*yd(i))/&
      &      1920.
          a0(i) = (9.*(y(i+2)+y(i-2))-116.*(y(i+1)+y(i-1))+2134.*y(i))/1920.
          a1d(i) = (34.*(yd(i+1)-yd(i-1))-5.*(yd(i+2)-yd(i-2)))/384.
          a1(i) = (-(5.*(y(i+2)-y(i-2)))+34.*(y(i+1)-y(i-1)))/384.
          a2d(i) = (12.*(yd(i+1)+yd(i-1))-yd(i+2)-22.*yd(i)-yd(i-2))/384.
          a2(i) = (-y(i+2)+12.*(y(i+1)+y(i-1))-22.*y(i)-y(i-2))/384.
          a3d(i) = (yd(i+2)-2.*(yd(i+1)-yd(i-1))-yd(i-2))/768.
          a3(i) = (y(i+2)-2.*(y(i+1)-y(i-1))-y(i-2))/768.
          a4d(i) = (yd(i+2)-4.*(yd(i+1)+yd(i-1))+6.*yd(i)+yd(i-2))/3840.
          a4(i) = (y(i+2)-4.*(y(i+1)+y(i-1))+6.*y(i)+y(i-2))/3840.
        END DO
        a0d(n-1) = (26.*yd(n-1)-yd(n)-yd(n-2))/24.
        a0(n-1) = (26.*y(n-1)-y(n)-y(n-2))/24.
        a1d(n-1) = (yd(n)-yd(n-2))/16.
        a1(n-1) = (y(n)-y(n-2))/16.
        a2d(n-1) = (yd(n)+yd(n-2)-2.*yd(n-1))/48.
        a2(n-1) = (y(n)+y(n-2)-2.*y(n-1))/48.
        a3d(n-1) = 0.0
        a3(n-1) = 0.
        a4d(n-1) = 0.0
        a4(n-1) = 0.
        wd(1) = 0.0
        w(1) = 1.
        wd(n) = 0.0
        w(n) = 1.
        fmd(n) = 0.0
        fm(n) = 0.
        IF (0. .GT. c(n-1)) THEN
          amin11 = c(n-1)
        ELSE
          amin11 = 0.
        END IF
        cl = -amin11
        IF (y(n) .GT. cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5)) THEN
          fmd(n-1) = cl*(yd(n)-(1.-cl)*0.5*(yd(n)-yd(n-1)))
          fm(n-1) = cl*(y(n)-(1.-cl)*(y(n)-y(n-1))*0.5)
        ELSE
          fmd(n-1) = yd(n)
          fm(n-1) = y(n)
        END IF
        clm = cl
        DO i=n-1,2,-1
          cl = clm
          IF (0. .GT. c(i-1)) THEN
            amin12 = c(i-1)
          ELSE
            amin12 = 0.
          END IF
          clm = -amin12
          x1 = 1. - 2.*cl
          x2 = x1*x1
          x3 = x1*x2
          IF (y(i) .GT. y(i+1)) THEN
            ymind = yd(i+1)
            ymin = y(i+1)
          ELSE
            ymind = yd(i)
            ymin = y(i)
          END IF
          IF (y(i) .LT. y(i+1)) THEN
            ymaxd = yd(i+1)
            ymax = y(i+1)
          ELSE
            ymaxd = yd(i)
            ymax = y(i)
          END IF
          IF (0. .LT. a0(i)*cl - a1(i)*(1.-x2) + a2(i)*(1.-x3) - a3(i)*(1.-x1*&
      &        x3) + a4(i)*(1.-x2*x3)) THEN
            fmimd = cl*a0d(i) - (1.-x2)*a1d(i) + (1.-x3)*a2d(i) - (1.-x1*x3)*&
      &        a3d(i) + (1.-x2*x3)*a4d(i)
            fmim = a0(i)*cl - a1(i)*(1.-x2) + a2(i)*(1.-x3) - a3(i)*(1.-x1*x3)&
      &        + a4(i)*(1.-x2*x3)
          ELSE
            fmim = 0.
            fmimd = 0.0
          END IF
          IF (fmim .GT. y(i) - ymin + fm(i)) THEN
            fmimd = yd(i) - ymind + fmd(i)
            fmim = y(i) - ymin + fm(i)
          ELSE
            fmim = fmim
          END IF
          IF (fmim .LT. y(i) - ymax + fm(i)) THEN
            fmimd = yd(i) - ymaxd + fmd(i)
            fmim = y(i) - ymax + fm(i)
          ELSE
            fmim = fmim
          END IF
          IF (0. .LT. fmim - (cl-clm)*y(i)) THEN
            fmd(i-1) = fmimd - (cl-clm)*yd(i)
            fm(i-1) = fmim - (cl-clm)*y(i)
          ELSE
            fmd(i-1) = 0.0
            fm(i-1) = 0.
          END IF
        END DO
        IF (0. .LT. c(1)) THEN
          cr = c(1)
        ELSE
          cr = 0.
        END IF
        y1d = cr*(yd(1)+(1.-cr)*0.5*(yd(2)-yd(1)))
        y1 = cr*(y(1)+(1.-cr)*(y(2)-y(1))*0.5)
        IF (y(1) .GT. y1) THEN
          fpd(1) = y1d
          fp(1) = y1
        ELSE
          fpd(1) = yd(1)
          fp(1) = y(1)
        END IF
        crp = cr
        wd(1:n) = 0.0
        DO i=2,n-1
          cr = crp
          IF (0. .LT. c(i)) THEN
            crp = c(i)
          ELSE
            crp = 0.
          END IF
          x1 = 1. - 2.*cr
          x2 = x1*x1
          x3 = x1*x2
          IF (y(i-1) .GT. y(i)) THEN
            ymind = yd(i)
            ymin = y(i)
          ELSE
            ymind = yd(i-1)
            ymin = y(i-1)
          END IF
          IF (y(i-1) .LT. y(i)) THEN
            ymaxd = yd(i)
            ymax = y(i)
          ELSE
            ymaxd = yd(i-1)
            ymax = y(i-1)
          END IF
          y2d = cr*a0d(i) + (1.-x2)*a1d(i) + (1.-x3)*a2d(i) + (1.-x1*x3)*a3d(i&
      &      ) + (1.-x2*x3)*a4d(i)
          y2 = a0(i)*cr + a1(i)*(1.-x2) + a2(i)*(1.-x3) + a3(i)*(1.-x1*x3) + &
      &      a4(i)*(1.-x2*x3)
          IF (0. .LT. y2) THEN
            fpid = y2d
            fpi = y2
          ELSE
            fpi = 0.
            fpid = 0.0
          END IF
          IF (fpi .GT. y(i) - ymin + fp(i-1)) THEN
            fpid = yd(i) - ymind + fpd(i-1)
            fpi = y(i) - ymin + fp(i-1)
          ELSE
            fpi = fpi
          END IF
          IF (fpi .LT. y(i) - ymax + fp(i-1)) THEN
            fpid = yd(i) - ymaxd + fpd(i-1)
            fpi = y(i) - ymax + fp(i-1)
          ELSE
            fpi = fpi
          END IF
          IF (0. .LT. fpi - (cr-crp)*y(i)) THEN
            fpd(i) = fpid - (cr-crp)*yd(i)
            fp(i) = fpi - (cr-crp)*y(i)
          ELSE
            fpd(i) = 0.0
            fp(i) = 0.
          END IF
          IF (fm(i-1) + fp(i) + 1.e-15 .LT. y(i)) THEN
            amax11d = yd(i)
            amax11 = y(i)
          ELSE
            amax11d = fmd(i-1) + fpd(i)
            amax11 = fm(i-1) + fp(i) + 1.e-15
          END IF
          wd(i) = (yd(i)*amax11-y(i)*amax11d)/amax11**2
          w(i) = y(i)/amax11
        END DO
        DO i=2,n-1
          yd(i) = yd(i) - (fmd(i-1)+fpd(i))*w(i) - (fm(i-1)+fp(i))*wd(i) + fmd&
      &      (i)*w(i+1) + fm(i)*wd(i+1) + fpd(i-1)*w(i-1) + fp(i-1)*wd(i-1)
          y(i) = y(i) - (fm(i-1)+fp(i))*w(i) + fm(i)*w(i+1) + fp(i-1)*w(i-1)

           WRITE(0,945) i, &
               y(i-1),  yd(i-1),  &
               y(i),    yd(i),    &
               y(i+1),  yd(i+1),  &
               fm(i-1), fmd(i-1), &
               fm(i),   fmd(i-1), &
               fm(i+1), fmd(i+1), &
               fp(i-1), fpd(i-1), &
               fp(i),   fpd(i),   &
               fp(i+1), fpd(i+1), &
               w(i-1),  wd(i-1),  &
               w(i),    wd(i),    &
               w(i+1),  wd(i+1); 
        END DO

        RETURN

945    FORMAT(//1x, 'i=', I5 // &
            4x, 'y(i-1)  = ', f5.3, 10x, 'yd(i-1)  = ',  f5.3    /  &
            4x, 'y(i)    = ', f5.3, 10x, 'yd(i)    = ',  f5.3    /  &
            4x, 'y(i+1)  = ', f5.3, 10x, 'yd(i+1)  = ',  f5.3    // &
            4x, 'fm(i-1) = ', f5.3, 10x, 'fmd(i-1) = ', f5.3    /  &
            4x, 'fm(i)   = ', f5.3, 10x, 'fmd(i)   = ', f5.3    /  &
            4x, 'fm(i+1) = ', f5.3, 10x, 'fmd(i+1) = ', f5.3    // &
            4x, 'fp(i-1) = ', f5.3, 10x, 'fpd(i-1) = ', f5.3    /  &
            4x, 'fp(i)   = ', f5.3, 10x, 'fpd(i)   = ', f5.3    /  &
            4x, 'fp(i+1) = ', f5.3, 10x, 'fpd(i+1) = ', f5.3    // &
            4x, 'w(i-1)  = ', f5.3, 10x, 'wd(i-1)  = ', f5.3    /  &
            4x, 'w(i)    = ', f5.3, 10x, 'wd(i)    = ', f5.3    /  &
            4x, 'w(i+1)  = ', f5.3, 10x, 'wd(i+1)  = ', f5.3    //)

      END SUBROUTINE ADV4M_D
